#include "Astronomy/CDynamicalTime.hpp"
#include <array>

// TD - UT1 计算表;
constexpr std::array<double, 112> dt_at =
{{
-4000,108371.7,-13036.80,392.000, 0.0000,
 -500, 17201.0,  -627.82, 16.170,-0.3413,
 -150, 12200.6,  -346.41,  5.403,-0.1593,
  150,  9113.8,  -328.13, -1.647, 0.0377,
  500,  5707.5,  -391.41,  0.915, 0.3145,
  900,  2203.4,  -283.45, 13.034,-0.1778,
 1300,   490.1,   -57.35,  2.085,-0.0072,
 1600,   120.0,    -9.81, -1.532, 0.1403,
 1700,    10.2,    -0.91,  0.510,-0.0370,
 1800,    13.4,    -0.72,  0.202,-0.0193,
 1830,     7.8,    -1.81,  0.416,-0.0247,
 1860,     8.3,    -0.13, -0.406, 0.0292,
 1880,    -5.4,     0.32, -0.183, 0.0173,
 1900,    -2.3,     2.06,  0.169,-0.0135,
 1920,    21.2,     1.69, -0.304, 0.0167,
 1940,    24.2,     1.22, -0.064, 0.0031,
 1960,    33.2,     0.51,  0.231,-0.0109,
 1980,    51.0,     1.29, -0.026, 0.0032,
 2000,    63.87,    0.1,   0,     0,
 2005,    64.7,     0.21,  0,     0,
 2012,    66.8,     0.22,  0,     0,
 2018,    69.0,     0.36,  0,     0,
 2028,    72.6
}};


double CDynamicalTime::Extrapolated(double year, double leapSecond) noexcept
{
    double dy = (year - 1820) / 100.0;
    return  -20.0 + leapSecond * dy * dy;
}

double CDynamicalTime::DeltaT(double y) noexcept
{
    ////表中最后一年;
    double y0 = dt_at[dt_at.size() - 2];
    //表中最后一年的deltatT;
    double t0 = dt_at[dt_at.size() - 1];
    // 超出部分用公式推算;
    if(y >= y0)
    {
        //leapSecond是y1年之后的加速度估计。瑞士星历表leapSecond=31,NASA网站leapSecond=32,skmap的leapSecond=29;
        double leapSecond = 31;
        if(y > y0 + 100.0)
            return Extrapolated(y, leapSecond);

        // 俩次推算,用来提高经度;
        //二次曲线外推;
        double v = Extrapolated(y, leapSecond);
        //ye年的二次外推与te的差;
        double dv = Extrapolated(y0, leapSecond) - t0;
        return v - dv * (y0 + 100.0 - y)/100.0;
    }
    //查表;
    int i = 0;
    for (; i < dt_at.size(); i+=5) 
    {
        // 判断年所在的区间;
        if(y < dt_at[i + 5])
            break;
    }
    //三次插值， 保证精确性;
    double t1 = (y - dt_at[i]) / (dt_at[i + 5] - dt_at[i]) * 10.0;
    double t2 = t1 * t1;
    double t3 = t2 * t1;

    return dt_at[i + 1] + dt_at[i + 2] * t1 + dt_at[i + 3] * t2 + dt_at[i + 4] * t3;
}
double CDynamicalTime::TDUT(double T) noexcept
{
    return DeltaT(T / 365.2425 + 2000.0)/86400.0;
}
